---
title: "sa_na_enmbg"
author: "nicheerT"
date: "2020/9/17"
output: html_document
---
### 准备数据：
```{r warning=FALSE}
library(ecospat)
library(ade4)
library(raster)
library(rworldmap)
Sys.time()
# C:\Users\Administrator\Desktop\envsV4tif
## 加载物种分布数据：
na_occ <- read.csv("D:/XH/xh/na.csv")[,2:3]
sa_occ <- read.csv("D:/XH/xh/sa.csv")[,2:3]
## 加载环境背景数据：
# 加载bg采样点：
#### 方法1： bg_enm
na_e_bg <- read.csv("D:/XH/fifth_bg2/bg_enm_buffer/na_enm_bg.csv")[,2:3]
sa_e_bg <- read.csv("D:/XH/fifth_bg2/bg_enm_buffer/sa_enm_bg.csv")[,2:3]
## 加载栅格背景数据：
env <- stack(list.files("D:/XH/third_env/envsV5tif",pattern = "tif",full.names = T))
构建数据组：
```{r warning=FALSE} env_na_occ <- data.frame(extract(env,na_occ)) env_sa_occ <- data.frame(extract(env,sa_occ)) names(env_na_occ)
方法1：bg_enm：
env_na_e <- data.frame(extract(env,na_e_bg)) env_sa_e <- data.frame(extract(env,sa_e_bg))
构建nat和inv：
构建nat_sa：
spocc <- c(rep(1,nrow(sa_occ)),rep(0,nrow(env_sa_e ))) envs_enm <- rbind(env_sa_occ,env_sa_e) names(sa_e_bg) <- names(sa_occ) sa_bg_occ <- rbind(sa_occ,sa_e_bg) nat_sa <- na.exclude(cbind(sa_bg_occ,envs_enm ,spocc))
构建inv_na:
spocc3 <- c(rep(1,nrow(na_occ)),rep(0,nrow(env_na_e ))) envs_enm3 <- rbind(env_na_occ,env_na_e) names(na_e_bg) <- names(na_occ) na_bg_occ3 <- rbind(na_occ,na_e_bg) inv_na <- na.exclude(cbind(na_bg_occ3,envs_enm3 ,spocc3))
names(inv_na) <- names(nat_sa)
```
sa_na_enmbg_PCA
```{r warning=FALSE}
生态位评估：第一阶段使用方法1：bg_enm：
pca.env <-dudi.pca(rbind(nat_sa,inv_na)[3:13], center = T, scale = T, scannf = F, nf = 2)
ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
PCA结果可视化：
library(ade4) library(factoextra) library(FactoMineR)
查看对应轴的贡献：
fviz_cos2(pca.env, choice = "var", axes = 1) fviz_cos2(pca.env, choice = "var", axes = 2)
重新构建矩阵可视化分布点：
saenm <- cbind(nat_sa[1:nrow(na.exclude(env_sa_occ)),3:13],"sa_occ") naenm <- cbind(inv_na[1:nrow(na.exclude(env_na_occ)),3:13],"na_occ") sabgenm <- cbind(nat_sa[-c(1:nrow(na.exclude(env_sa_occ))),3:13],"sa_bg") nabgenm <- cbind(inv_na[-c(1:nrow(na.exclude(env_na_occ))),3:13],"na_bg") names(saenm) <- names(naenm) <- names(sabgenm) <- names(nabgenm) sdm3 <- rbind(sabgenm,nabgenm,saenm,naenm) colnames(sdm3)[12] <- "group" sdm3$group <- factor(sdm3$group,levels=c("sa_bg","na_bg","sa_occ","na_occ") ) pca3 <- PCA(sdm3[,1:11],graph = FALSE) fviz_pca_biplot(pca3, col.ind = sdm3$group, palette = c("pink","lightskyblue","red","blue"), addEllipses = FALSE, label = "var", col.var = "black", repel = TRUE, legend.title = "Species",title ="SA_na_PCA") ```
COMP_NICEH
```{R warning=FALSE}
predict the scores on the PCA axes
scores.globclim <- pca.env$li
PCA scores for the species native distribution
scores.sp.nat <- suprow(pca.env,nat_sa[which(nat_sa[,14]==1),3:13])$li
PCA scores for the species invasive distribution
scores.sp.inv <- suprow(pca.env,inv_na[which(inv_na[,14]==1),3:13])$li
PCA scores for the whole native study area
scores.clim.nat <- suprow(pca.env,nat_sa[,3:13])$li
PCA scores for the whole invaded study area
scores.clim.inv <- suprow(pca.env,inv_na[,3:13])$li
calculation of occurence density
z1<- ecospat.grid.clim.dyn(glob=scores.globclim, glob1=scores.clim.nat, sp=scores.sp.nat, R=100,th.sp=0.01) z2 <- ecospat.grid.clim.dyn(glob=scores.globclim, glob1=scores.clim.inv, sp=scores.sp.inv, R=100,th.sp=0.01)
overlap corrected by availabilty of background conditions
ecospat.niche.overlap(z1,z2,cor=T)
uncorrected overlap
ecospat.niche.overlap(z1,z2,cor=F)
test of niche equivalency and similarity
test of niche equivalency
Niche equivalency test H1: Is the overlap between
the native and invaded niche higher than two random niches?
查看cpus
library(parallel) source("C:/Users/admin/Desktop/ecospat.niche.equivalency.test.R") source("C:/Users/admin/Desktop/ecospat.niche.similarity.test.R") eq.test1 <- ecospat.niche.equivalency.test2(z1, z2, rep=100, alternative = "lower",ncores = 10) eq.test1 eq.test <- ecospat.niche.equivalency.test2(z1, z2, rep=100, alternative = "greater",ncores = 10) eq.test$p.D eq.test$p.I
ecospat.plot.overlap.test(eq.test1, "D", "Equivalency")
test of niche similarity
sim1<-ecospat.niche.similarity.test2(z1,z2,rep=100,alternative = "lower",rand.type = 1,ncores = 10) sim1 # sim2<-ecospat.niche.similarity.test2(z1,z2,rep=100,alternative = "greater",rand.type = 2,ncores = 10) sim2$p.D sim2$p.I
#niche randomly shifted only in invaded area
ecospat.plot.overlap.test(sim2, "D", "Similarity")
绘制可视化重叠：
ecospat.plot.niche.dyn(z1, z2, quant=0.25, interest=2, title= "SA vs AS Niche Overlap", name.axis1="PC1(tem)",name.axis2="PC2(pre)") ecospat.shift.centroids(scores.sp.nat, scores.sp.inv, scores.clim.nat, scores.clim.inv)
观测比例：
sa_na_index <- ecospat.niche.dyn.index (z1, z2,intersection=0.01) sa_na_index$dynamic.index.w
sa_na_index <- ecospat.niche.dyn.index (z1, z2,intersection=0.05) sa_na_index$dynamic.index.w
sa_na_index <- ecospat.niche.dyn.index (z1, z2,intersection=0) sa_na_index$dynamic.index.w
sa_na_index <- ecospat.niche.dyn.index (z1, z2,intersection=NA) sa_na_index$dynamic.index.w ```
```